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Assessing the robustness of spatial pattern 
sequences in a dryland vegetation model 

Kama Gowda* Yuxin Chen* Sarah lamsj Mary Silber^ 


A particular sequence of patterns, “gaps —)■ labyrinth —)■ spots,” occurs with de¬ 
creasing precipitation in previously reported numerical simulations of PDE dryland 
vegetation models. These observations have led to the suggestion that this sequence 
of patterns can serve as an early indicator of desertification in some ecosystems. Since 
parameter values can take on a range of plausible values in the vegetation models, it 
is important to investigate whether the pattern sequence prediction is robust to varia¬ 
tion. For a particular model, we find that a quantity calculated via bifurcation-theoretic 
analysis appears to serve as a proxy for the pattern sequences that occur in numerical 
simulations across a range of parameter values. We find in further analysis that the 
quantity takes on values consistent with the standard sequence in an ecologically rel¬ 
evant limit of the model parameter values. This suggests that the standard sequence 
is a robust prediction of the model, and we conclude by proposing a methodology for 
assessing the robustness of the standard sequence in other models and formulations. 


1. Introduction 

Many studies of spatially periodic patterns in models of dryland vegetation focus on patterns 
as potential indicators of ecosystem vitality [1-13]. In particular, flat-terrain patterned states 
in several models have been shown in simulations to evolve through a sequence of morphologies, 
“gaps —)■ labyrinth —)■ spots” (Figure 1), as ecosystem aridity increases [1, 4, 14, 15]. This sequence 
precedes the collapse of vegetation in the models, which has led to the suggestion that real ecosys¬ 
tems may evolve through this predictable sequence of patterns en route to desertification [16, 17]. 
In this way, vegetation patterns may serve as early-warning signs of critical ecosystem transitions. 

The pattern sequence prediction emerges from a modelling framework comprising a number of 
different ecological hypotheses, functional formulations, and restrictions on plausible parameter sets. 
It is therefore important to investigate whether the prediction is robust within this framework. It 
is encouraging that multiple models can produce the pattern sequence “gaps —)■ labyrinth —)■ spots” 
(the standard sequence), but little work has been done to assess whether the sequence occurs in any 
model under variations of the parameter values, which may vary significantly between ecosystems. 
In this paper, we ask whether easily calculable quantities, derived from bifurcation theoretic analysis 
of pattern-forming instabilities, can be used to predict where the standard sequence occurs in a 
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Figure 1: Example of the standard “gaps —> labyrinth —>■ spots” sequence in the vegetation model by Rietkerk 
et al. Qualitatively different patterns occur at successively smaller values of a precipitation 
parameter. Darker shading denotes higher levels of vegetation biomass. 


model’s parameter space. To do this, we analyse and simulate the widely studied model by Rietkerk 
et al. [14] across a broad range of parameter values. 

The standard sequence, “gaps —)■ labyrinth —)■ spots”, is observed in a suite of partial differential 
equation (PDF) models that describe the community-scale dynamics of dryland vegetation over 
flat terrain [1, 4, 14, 15]. These models all involve a Turing instability [18] in the formation of 
patterns from uniform vegetation. Observations of the standard sequence as a model prediction 
come primarily from numerical simulations [1, 4, 14], which show gap, labyrinth, and spot patterns 
occurring at successively lower values of a precipitation parameter (Figure 1). A study by LeJeune et 
al. [15] used bifurcation analysis to demonstrate analytically that this sequence occurs in a tractable 
1-field model for a particular parameter set. The apparent agreement between these observations, 
which come from different model formulations, provides some support for the standard sequence as 
a robust prediction of this suite of models. 

Empirical support for this sequence comes chiefly from two studies of remotely-sensed imagery. 
A 2006 study by Barbier et al. [19] used imagery over southwest Niger to demonstrate that gap 
patterns emerged from uniform vegetation over a period of time coinciding with a prolonged Sahe¬ 
lian drought. This result is consistent with model observations of gap patterns occurring near the 
onset of pattern formation. A 2011 study by Deblauwe et al. [20] classified pattern morphologies 
in imagery over Sudan and found that different morphologies vary over spatial precipitation gra¬ 
dients in accordance with the standard sequence prediction. Gaps tended to occur in areas with 
relatively high mean annual precipitation, spots occurred in areas with relatively low precipita¬ 
tion, and labyrinths occurred in between. Pattern dynamics were also assessed using three sets 
of images taken over a 35-year span. Gaps in some areas were shown to transition to labyrinths 
over a period of time again coinciding with a sustained regional drought. Labyrinths transitioned 
to spots in different areas over the same period of time. Though neither of these studies show 
the standard sequence preceding the collapse of vegetation, they demonstrate consistency between 
model predictions and empirical observations. 

Since these empirical studies focus on two Sub-Saharan regions, they make statements about 
particular ecosystems, which are subject to particular environmental and ecological parameters. It 
is unknown how other periodically-patterned flat-terrain ecosystems behave in response to changes 
in aridity. It is also unknown whether the suite of dryland vegetation models predicts the standard 
sequence to occur under most circumstances, since the numerical simulations underlying the stan¬ 
dard sequence model observations were generated using a small number of parameter sets. The 
parameter sets appropriate for different ecosystems can differ significantly. For instance, the model 
for banded vegetation patterns on slopes by Klausmeier [21] distinguishes between shrublands and 
grasslands, using plant mortality parameters that differ by an order of magnitude for these differ- 
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ent plant types^. Little work has been done to systematically assess the sensitivity of the standard 
sequence prediction to model parameters. Such an assessment is important to understanding under 
what circumstances the sequence is predicted to be an early-warning sign. 

Gowda et al. [5] demonstrate that, in certain limits of model parameter values, pattern sequences 
can be studied analytically using bifurcation theory. They compute the coefficients of equations 
which describe the amplitudes of Fourier modes on a 2D hexagonal lattice near a pattern-forming 
instability. These equations can provide a reduced description of dynamics near the instability. In 
limits of the model parameter values where transitions between patterns all occur in a regime of weak 
nonlinearity, their analysis shows that the coefficient of the quadratic-order term in these equations 
affects the sequence that occurs. If this coefficient changes its sign from negative to positive as a 
precipitation parameter is decreased in value, an analog of the standard sequence occurs in certain 
cases. Otherwise, alternative sequences, such as one consisting only of spot patterns, can occur. 
Based on a preliminary numerical investigation of the model by von Hardenberg et al. [1], Gowda et 
al. [5] speculate that this coefficient also encodes information about pattern sequences that occur 
in more fully nonlinear cases. 

Here, we ask whether the standard sequence can be identified in the model by Rietkerk et al. [14] 
using the sign of the quadratic coefficient of the hexagonal lattice amplitude equations. Specifically, 
we investigate whether the standard sequence occurs at parameter values where the quadratic 
coefficient changes its sign from negative to positive as the precipitation parameter decreases. We 
calculate the amplitude equation coefficients across a range of parameter values, and evaluate 
these coefficients at two values of the precipitation parameter which correspond to pattern-forming 
bifurcation points. We conduct numerical simulations of the model over the same parameter range 
in order to identify the pattern sequences that occur. We compare the results of the analysis with 
those of the simulations in order to address whether the standard sequence in this model is signalled 
by the quadratic coefficient changing signs as the precipitation parameter decreases. Of particular 
interest is whether this is true in regimes where the weakly nonlinear analysis provides no direct 
information about pattern stability. 

If coefficients from bifurcation analysis can serve as proxies for the standard sequence, then this 
would allow for more efficient exploration of a model’s parameter space than by direct numerical 
simulation. This is important for assessing the robustness of the standard sequence in 3-field 
models, such as the models by Rietkerk et al. [14] and Gilad et al. [11, 22], which depend on a 
large number of non-dimensional parameters. In this paper, we find that a natural limit of the 
Rietkerk et al. model parameter values permits analytical insight into the model terms which set 
the sign of the quadratic coefficient. This limit corresponds to a surface water diffusion rate which 
is much larger than the rate of surface water infiltration as well as the rate of biomass dispersal. 
Studying the model in this limit allows us to comment on the implications of our numerical study 
on an ecologically relevant portion of the model parameter space. 

This paper is organised as follows. Section 2 gives background information about the model by 
Rietkerk et al. and the bifurcation-theoretic analysis used in our study of patterned states. Section 3 
describes the methods used in our study, specifying the parameter spaces explored and summarising 
the analytic and numerical methods used to explore these parameter spaces. Section 4 describes the 
results of the analytic and numerical parameter study, as well as the results of additional analysis 
performed in a distinguished parameter limit (described in the appendices). Section 5 discusses the 
implications of our study. 


^Klausmeier [21] uses the mortality rate Mtree = 0.18 year ^ for shrublands and Mgrass = 1.8 year ^ for grasslands. 
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2. Background 


2.1. Model by Rietkerk et al. (2002) 


We study the PDE vegetation model by Rietkerk et al. [14] (R02), which consists of three fields: 
surface water /i, soil water w, and plant biomass n. Using the non-dimensional form given by 
Zelnik et al. [11], the model is written as 


where 


dth = p — I{n)h +DhV'^h, 

precip. infil. diffusion 

dtw = — + I{n)h — '^G{w)n + 

evap. transp. diffusion 

dfU = — fill +G{w)n+ , 

mort. growth dispersal 


I{n) = a —and G{w) = 


n + 1 


w + 1 


( 1 ) 


In this model, precipitation is a constant input to the surface water field. Surface water infiltrates 
and becomes soil water. The infiltration rate (i.e. the conversion rate of surface water to soil water) 
increases in the presence of biomass via the function I{n) to model the increased permeability of the 
soil due to plant roots, /(n) saturates to a as n —)■ oo. Water leaves the soil via evaporation, and is 
also transpired by plants. The growth rate of biomass is directly proportional to the transpiration 
rate, and increases with the availability of soil water via the saturating function G{w). Together, 
these terms make a positive feedback between infiltration and biomass growth: biomass growth 
increases with soil water, soil water increases with infiltration, and the infiltration rate increases 
with biomass. 

This model includes surface and soil water diffusion terms. Plant dispersal, which encompasses 
seed dispersal and clonal growth, is also modelled using a diffusion term. The diffusion terms are 
in two spatial dimensions (2D), i.e. Jdx^ + d‘^jd'ip'. Surface water diffusion is typically 

assumed to occur much more rapidly than soil water diffusion, so ^ Among three-field 
PDE vegetation models, soil water diffusion and plant dispersal have been modelled as occurring 
on either similar [14] or different [22] length and time scales with > 1. An advection term 
present in the original form of R02 is neglected here, because the focus of our investigation is on 
flat-terrain patterns. The dynamics of water on a slope modelled via advection break the symmetry 
that causes 2D patterns such as gaps or spots at pattern onset. 

In general, the form of the growth term varies between models [21, 22], and it determines the 
number of uniform steady state solutions that occur for a given system. Eor R02, the rate of biomass 
growth depends linearly on the amount of biomass. The growth rate is also a saturating function 
of soil water, so that it is linear in the amount of soil water for small values of this variable, and 
constant for large values. This growth rate permits two spatially uniform steady state solutions, 
which satisfy the equations 


0 = » — I(n)h = p — a — 

^ ^ n + l 

Ti ~\~ f VO 

0 = — vw -|- I{n)h — ^G{w)n = —uw + a - — 7- zru. 


0={-p + G{w)) n = (^-p + • 


n -|- 1 


rc -|- 1 
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Precipitation {p) 


Figure 2: Schematic diagram depicting the uniform steady state solutions of R 02 ( 1 ), with insets showing 
examples of patterned states occurring at different values of p. The uniform desert state is 
stable on the interval p € [ 0 ,po)j s-nd the uniform vegetated state is stable to spatially uniform 
perturbations for p > po. The vegetated state is unstable to spatially periodic perturbations at 
a range of wavelengths on the interval p € {pi,pu)- We refer to the endpoints pi and Pu as the 
lower and upper Turing points, respectively. 


One solution, {h,w,n) = {p/fa,p/ 1 ', 0 ), represents a zero-biomass desert state. The other solution 
{h,w,n) = {hQ,wo,nQ) represents a vegetated state with nonzero biomass, 

ho = 77 —r, Wo = - -, no = — [p- - - , ( 2 ) 

/(no) I- p. 'yp \ 1- pj 

for which no > 0 when p > pv/{l — p) = po- We note that the soil water wq in the vegetated state 
depends only on the mortality parameter p, and not on the precipitation level p. A diagram of 
uniform steady state biomass as a function of the precipitation parameter is shown in Figure 2. The 
desert state is stable to spatially uniform perturbations at low values of precipitation (0 < p < po)) 
while the vegetated state is stable to such perturbations at higher values of precipitation (p > po)- 
These steady states exchange stability in a transcritical bifurcation at p = po- 

For a range of precipitation values, indicated by the interval p G {p£,Pu), tbe vegetated state 
may be unstable to spatially periodic perturbations at a range of wavelengths. When this is the 
case, pattern-forming instabilities occur at pi and pu via the Turing mechanism [18], and we refer 
to these points as the lower and upper Turing points respectively. These instabilities result in the 
formation of spatial patterns, which are eventually stabilised by nonlinearities in the system. 

As models of dryland vegetation-water interactions, R02 and related models are conceptual in 
the sense that they are derived using simple mathematical assumptions for the forms of rate and 
transport terms. The predictions of the model are not considered robust in regimes where the 
results of calculations are sensitive to a modest variation in the model parameters or in the exact 
form of the feedbacks used to model interactions. We explore R02 across a range of parameter 
values to assess the robustness of model behaviours to parameters. 

We consider variations of the infiltration rate parameters, / and a. As biomass increases, the 
infiltration rate approaches a. For low values of biomass, infiltration occurs at the lower rate af, 
with 0 < / < 1. The parameter / controls the strength of the infiltration feedback, where / = 0 
corresponds to maximal feedback and / = 1 corresponds to no feedback. Some soil types have high 
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infiltration rates and are unaffected by the presence of vegetation, while others have low infiltration 
rates that are more readily influenced. Because of this, it is appropriate to consider a range of 
/ and a values in a parameter exploration of R02. We also consider variations of the surface 
water diffusion rate, Dh. Turing instabilities in R02 occur when plant dispersal and surface water 
diffusion occur on different scales. The ratio between surface diffusion and plant dispersal rates, 
Dh, is central to pattern formation in this model, and thus it is varied over a wide range in the 
investigations in this paper. 

2.2. Analysis of patterns and transitions 

Pattern formation in R02 and related PDE models is a nonlinear phenomenon and can be studied 
analytically via bifurcation theory only in certain limits [23, 24]. One such limit occurs in the 
vicinity of a Turing point [18, 24]. At a Turing point, a uniform steady state has a zero linear growth 
rate (is neutrally stable) for Fourier mode perturbations of a particular wavelength. We refer to 
this wavelength as the critical wavelength Ac, corresponding to a critical wave number Qc = 27r/Ac. 
The uniform vegetated state of the R02 model can be unstable to Fourier mode perturbations over 
an interval delimited by two Turing points, which we refer to as the lower and upper Turing points, 
P£ and Pu respectively. In this investigation, we analyse patterned states via bifurcation theory near 
these two points by considering only Fourier modes of a solution on a 2D hexagonal lattice [23, 24]. 
The lattice is constructed so that Fourier modes associated with the critical wavelength grow in 
a neighbourhood of the Turing point, while all others are linearly damped [25]. A system of 
six ordinary differential equations (ODFs) for the Fourier mode amplitudes and their complex 
conjugates can be written to describe the local pattern-forming dynamics of these critical modes. 
We study the steady states of these ODFs to gain insight into a system’s pattern-forming behaviour 
near the Turing points. 

The critical Fourier modes on a hexagonal lattice take the form 

+ C.C., (3) 

where zi, Z 2 , and zs are time-dependent complex amplitudes and c.c. denotes the complex conjugates 
of these modes. The vectors qi, q 2 , and qs are wave vectors such that 

qi = qc(l, 0), q 2 = gc(-l/2, \/3/2), qa = -(qi -F q 2 ), 

where Qc is the critical wave number. By assuming that pattern dynamics near a Turing point 
are well-approximated by the modes (3) as a small perturbation to the local uniform steady state, 
truncated equations for the amplitude dynamics are found to be [23, 26] 


Zi 

= mzi -1- az 2 Z 3 - ( 

b\zi\^ -f c(l2;2l^ -f 

i^3p); 

) Zl, 


Z2 

= mz2 + aziZ3 - ( 

h\z2\^ + c(l2;il^ + 

i^3p); 

) Z2, 

(4) 


= mz3 -1- aziZ2 - ( 

b\z3\^ + c{\zi\‘^ + 

l^2p)' 

)zs. 



These equations are considered to be valid in a regime of weak nonlinearity in a neighbourhood of 
the Turing point. In this regime, the amplitudes are small so that all three terms can contribute to 
dynamics on the same order. The linear coefficient m is equivalent to the growth rate of the critical 
modes obtained by linearising the original system, so that it is small near the Turing point and 
exactly zero at that point. The quadratic coefficient a must be small for any steady state solutions 
to (4) to be stable [26], and must be on the order of the size of the amplitude for the quadratic term 
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Table 1: Equations for steady state solutions to (4), and distinct eigenvalues for these solutions, 
z = (zi, Z 2 , Z 3 ), and Xs, Xh > 0 . 



Branching equation 

Distinct eigenvalues 

Gaps {H~): 

z = -{Xh,Xh,Xh) 

2>axh, 2axh + 2(c - 


0 = m — axh — {b + 2c)x| 

—axh — 2(6 -|- 2c)x\, 0 

Stripes (S): 

z = (a:s,0,0) 

—26x^, —(c — b)Xg — axg, 


0 = m — hxl 

— (c — b)x‘l + axs, 0 

Spots {H^): 

Z = {xh,Xh,Xh) 

-3axh, -2axh + 2 (c - 6 )x|, 


0 = m + axh — (^ -|- 2 c)x| 

axh — 2(6 -|- 2c)x|, 0 



gaps {H ) stripes (5) spots {H^) 


Figure 3: Examples of hexagon {H and and stripe {S) patterns on a 2D hexagonal lattice. We idealise 
gaps as H~ patterns, labyrinths as S, and spots as . 


to balance with the linear and cubic terms. The cubic coefficients b and c are saturating terms, 
and also affect the stability of steady state solutions of (4). 

Hexagon and stripe steady state solutions of (4) are summarised in Table 1. Hexagon steady 
states correspond to equal, real-valued amplitudes, with gaps being negative and spots positive. 
Stripes correspond to the case of only one nonzero amplitude. The stability of these small-amplitude 
solutions to perturbations on the hexagonal lattice is dictated by the eigenvalues listed in Table 1. 
Of primary importance to this study is the result that a < 0 is a necessary condition for the stability 
of gaps, and a > 0 for spots. The quantities b, c — b and 6 -|- 2c are also related to stability and we 
refer to these in later interpretation of results. 

Gap, labyrinth, and spot-patterned states have previously been observed in R02 [14]. Examples 
of these states are shown in Figure 2. We idealise these states as hexagon and stripe steady states 
on the hexagonal lattice as in [5, 15] (Figure 3). In general, the lower and upper Turing points occur 
sufficiently far from one another that it is not possible to study transitions between these states via 
local bifurcation theory in scenarios of changing precipitation. However, Gowda et al. [5] consider 
the special limit in which two Turing points occur near to one another (near to a degenerate Turing 
point), resulting in pattern transitions that can be analysed via bifurcation theory. In a generic 
analysis, they find that an analog of the standard “gaps —)■ labyrinth —)■ spots” sequence is one of 
many possible sequences that can occur as precipitation decreases. They also find that the value 
of the quadratic coefficient a at the lower and upper Turing points may serve as a proxy for the 
apparent pattern sequence. Given appropriate conditions on the cubic coefficients b and c that 
allow the stability of hexagons and stripe solutions to (4), they find that if a < 0 at the upper 
Turing point, gaps appear at the start of the sequence. If a > 0 at the lower Turing point, then 
spots appear at the end of the sequence. If these two conditions occur together, then some analog 
of the standard sequence can occur. 
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Gowda et al. [5] speculated that the quadratic coefficient a can be used to obtain information 
about pattern sequences outside the limited setting of a degenerate bifurcation problem. They 
conducted a bifurcation analysis and small set of numerical simulations on the model by von 
Hardenberg et al. [1] (VHOl). VHOl is a two-field PDE vegetation model with cross-diffusion which 
features two Turing bifurcations on a uniform vegetated steady state (similar to R02). Varying a 
diffusion coefficient and leaving all other parameters of VHOl fixed, Gowda et al. [5] found that 
a < 0 at the upper Turing point and a > 0 at the lower Turing point in a regime where the 
diffusion coefficient is sufficiently large. Where stable solutions to the amplitude equations exist 
within this regime, the signs of the a coefficient at the Turing points allow pattern sequences 
to begin with stable gaps and to end with stable spots as precipitation decreases. In between 
these patterns, small-amplitude stripe patterns appear in both analysis and numerical simulations. 
Where the cubic coefficients b and c prevent the stability of small-amplitude patterns, disordered 
gaps, labyrinths and spots states were again observed with decreasing precipitation in a numerical 
simulation. This observation led Gowda et al. [5] to speculate that the standard sequence occurs 
in regimes of the model parameter space demarcated by the a coefficient taking a negative sign at 
the upper Turing point and a positive sign at the lower Turing point, regardless of whether the 
amplitude equations (4) predict stable patterned states in their small-amplitude regime of validity. 
Here, we examine whether this holds in the R02 model. 

3. Methods 

We ask whether the quadratic coefficient a of the amplitude equations (4) evaluated at the upper 
and lower Turing points signals the appearance of the standard sequence in the model by Rietkerk et 
al. [14] (R02). To answer this, we conducted bifurcation analyses and numerical simulations across 
two relevant parameter spaces of the model. To investigate which model terms and parameters 
influence the sign of the quadratic coefficient, we conducted additional analysis in an ecologically 
relevant distinguished limit of the model parameters. 

3.1. Model parameters 

We studied variations of the non-dimensional R02 parameters given by Zelnik et al. [11], which 
are based on the dimensioned parameter values estimated by Rietkerk et al. [14]. These parameter 
values are summarised in Table 2. The / parameter in R02 controls the strength of the infiltration 
feedback, and is bounded between 0 and 1. The default value given in [11] is / = 0.2. We used 
/ G [0.1, 0.9] in numerical simulations. The a parameter controls the rate of infiltration, and can 
plausibly take on a large range of values depending on the soil type that is modelled. The default 
value given in [11] is a = 0.4, and we considered a G [10“^, 10^]. The parameter is the ratio of 
the surface water diffusion rate and the biomass dispersal rate. Rietkerk et al. [14] use Dh = 10^ 
for R02 and Zelnik et al. [11] use = 10^ for the corresponding parameter value in a simplified 
version of the model by Gilad et al. [22]. We varied Dh G [10°'®, 10^]. To consider the dependence of 
results on co-variation of parameters, we studied the a-f and Df^-f parameter spaces. In additional 
analysis described in the appendices and summarised in Section 44.3, we considered an ecologically 
relevant limit where Df^ ^ a. Our parameter exploration was conducted primarily in this limit, 
since a was held fixed at 0.4 while Dh was varied, and Dh was held fixed at 10^ while a was varied. 



Table 2: Parameters given by Zelnik et al. [11] for the R02 model, and parameters varied in this 
study. 



Interpretation 

Value in [11] 

Constraints 

Variation studied 

h 

mortality rate 

0.5 

0 < /i < 1 

— 

a 

infiltration rate 

0.4 

a > 0 

10-1 - 10^ 

f 

infiltration feedback strength 

0.2 

0 </<l 

0.1 - 0.9 

V 

evaporation rate 

0.4 

n >Q 

— 

7 

transpiration rate 

0.1 

7 > 0 

— 


ratio of soil water diffusion rate to 
biomass diffusion rate 

1 

l<D^<Dh 

— 

Dh 

ratio of surface water diffusion rate to 
biomass diffusion rate 

10 ^ 

Dh 7 

100.5 _ iqA 


3.2. Amplitude equation calculations 

We computed coefficients of the amplitude equations (4) for R02 using the procedure outlined in 
Judd et al. [25], which takes a perturbative approach to obtaining these coefficients for a 2-field 
reaction diffusion system. These coefficients are written as expressions of the reaction functions and 
diffusion parameters. We adapted the Judd et al. [25] procedure for a three-field reaction-diffusion 
system to obtain expressions for the amplitude equation coefficients; the aspect of this procedure 
specific to the calculation of the quadratic coefficient a is described in Appendix A. The values 
of the amplitude equation coefficients a, b, and c are computed at the lower and upper Turing 
points via Mathematiea. These calculations are performed on a grid of points in the a-f and D^-f 
parameter spaces. The results of these calculations were verified at a few non-degenerate points 
using a centre manifold reduction approach (the general approach is described in [27]). A Nelder- 
Mead minimisation library in Mathematiea was used to find roots of the quantities a, b, c — b, and 
b + 2c, which are relevant to assessing the stability of solutions to (4). 

We also calculated the scaling behaviour of the quadratic coefficient a with respect to the quantity 
6 = I{no)/Dh, where I{n) = a(n-|-/)/(n-|-l), and no is the uniform vegetated biomass given by ( 2 ). 
These calculations are presented in the appendices. We found that the critical wave number and 
onset parameter value scale in distinct ways at each Turing point with respect to 6 when S is 
sufficiently small. We used these scaling relationships to calculate a closed-form expression for the 
leading order value of a at each Turing point. 

3.3. Numerical simulations 

We conducted numerical simulations to identify pattern sequences that occur as precipitation de¬ 
creases in the R02 model. We employed a numerical procedure which simulates an ecosystem 
undergoing a slow monotonic change in precipitation over time. These simulations were run using 
grids of parameter values covering regions of the a-f and D^-f parameter spaces. The procedure is 
outlined schematically in Figure 4, and is described in detail in the electronic supplementary mate¬ 
rial. For each set of parameter values, precipitation is incremented in small discrete steps, and the 
solution is allowed to reach a steady state between these steps in precipitation. The final state at 
the previous precipitation value is used as the initial condition for the new precipitation value. The 
precipitation increment step size was chosen based on the distance between the upper and lower 
Turing points, Pu—pe, so that approximately 30-100 end states were saved per simulation. Discrete 
steps were chosen instead of continuously varying precipitation to avoid transient effects, i.e. sim- 
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Figure 4: Diagram of numerical simulation procedure. Numerical simulations are run at discrete values of p 
marked by dots. The procedure is initialised with p just below the upper Turing point (green 
star) and run forward in time until a steady state is reached. Then p is stepped upward by a 
small increment and the simulation is once again allowed to reach steady state. This is repeated 
until patterns lose stability to a uniform vegetated state at p = (red point, right). Using 
the previous patterned steady state (blue square, right) as an initial condition, p is then stepped 
downward in the same way until patterns lose stability to a uniform state at p = p^_ (red point, 
left). Then p is incremented upward a final time, and the procedure terminates when patterns 
once again lose stability (red octagon). Note that pu+ and pi- do not necessarily coincide with the 
Turing points Pu and pi, because patterns may persist outside of the Turing instability interval. 


ulation results that are sensitive to the rate at which precipitation changes. Simulations were run 
using the exponential time differencing Runge-Kutta 4 (ETDRK4) scheme [28, 29] modified for 2D 
systems [30]. 

The procedure is constructed to run simulations over the interval of p where patterns are sta¬ 
ble and to identify any possible history-dependence (hysteresis) in the pattern sequences. This is 
accomplished by first incrementing p upward until patterns die out to yield a uniform state. We 
denote this point of pattern die-off as p = Pu+- Precipitation is then decremented, simulating a 
scenario in which an ecosystem slowly becomes more arid. This continues until patterns die out 
again, which yields another uniform state, which we denote p = p£_. Precipitation is incremented 
upward a final time to assess hysteresis in the pattern sequence (i.e. whether the sequence occurs 
differently when p is slowly increasing versus decreasing). The procedure terminates when patterns 
die out once more. An approximate interval for the stability of patterns is given by p G {p(,-,pu+), 
which contains the Turing instability interval p G {pi,Pu)- The Turing instability interval is deter¬ 
mined via a linear stability calculation and does not capture the nonlinear stabilisation of patterns. 
In cases where the amplitude equations (4) predict stable hexagons solutions near a Turing point, 
these solutions bifurcate in such a way as to be stable outside the Turing interval. When amplitude 
equation solutions branch away from the Turing interval (e.g. when stripes bifurcate subcritically), 
these solutions may also stabilise at large amplitude outside the Turing interval. 

4. Results 

4.1. Amplitude equation calculations 

We find that the a-f and Dh-f parameter spaces of the model by Rietkerk et al. [14] (R02) can be 
divided into regions where the amplitude equations (4) give different qualitative predictions. The 
results of the calculations at the upper Turing point are summarised in Figure 5, and the lower 
Turing point calculations are summarised in Figure 6. In the unlabelled white regions, no Turing 
points occur because the uniform vegetated steady state is stable to spatially periodic perturbations, 
and no calculations are performed. The black curves separating the white and shaded regions denote 
a degeneracy of the Turing points, where the upper and lower Turing points come together at a 
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single precipitation value. In the shaded regions below this curve, two Turing points occur on the 
vegetated state, and analysis is performed at each point. 

Each shaded region in Figures 5 and 6 is associated with a qualitatively distinct bifurcation 
diagram applicable to a neighbourhood of the Turing point. The qualitative aspects (e.g. stability, 
branching direction) of the bifurcation diagrams are determined by the signs of the quantities a, 6, 
c — b, and b + 2c. These quantities arise from the amplitude equation steady state eigenvalues and 
branching equations listed in Table 1. Notably, the sign of the a coefficient serves as a necessary 
condition for the stability of either small-amplitude gaps or spots solutions. A necessary condition 
for the stability of gaps is given by a < 0, and a necessary condition for the stability of spots is 
given by a > 0. 

For the analyses summarised in Figures 5 and 6, the regions are arrayed similarly in both a-f 
and Dh-f parameter spaces. For example, regions A-E in Figure 5 corresponding to the upper 
Turing point occur in the same order when varying parameters away from the Turing degeneracy 
curve (e.g. when increasing Dh compared to decreasing a, with / fixed). We find that this occurs 
because the Turing point calculation and the coefficients a, 6, and c depend on the quantity a/D^, 
and not a. and independently. In the appendices, we show how the Turing point calculation 
and the quadratic coefficient a depend on ajD^ via the quantity S oc a/D^. Though the amplitude 
equation coefficients do not depend on a and Dh independently, this does not translate to the 
invariance of full solutions to R02 for fixed ratios of a/Dh- This can be seen in weakly nonlinear 
solutions, which have linear eigenfunctions that depend on a and Dh independently. 

We first interpret the results of bifurcation analysis at the upper Turing point, which are sum¬ 
marised in Figure 5. We consider a scenario in which precipitation decreases slowly over time, so 
that the upper Turing point threshold is crossed from above. The sequence of pattern morphologies 
observed in such a scenario begins with patterns born near the upper Turing point. The regions in 
Figure 5 specify whether the amplitude equations predict a stable patterned state in some neigh¬ 
bourhood of the upper Turing point, and also the morphology of that state. A gap {H~) patterned 
state stable near the upper Turing point accords with the standard “gaps —)■ labyrinth —)■ spots” 
sequence prediction. 

In region A of Figure 5, the quantities a, c — b, b + 2c, and b are all positive, which allows a 
stable spot solution {H~^) to the amplitude equations in a neighbourhood of the upper Turing 
point. This analysis predicts that pattern sequences begin with spot patterns in region A of the 
parameter space, which is inconsistent with the standard sequence. The stripes solution (S) can 
also be stable in region A. It stabilises away from the Turing point, so that spots may transition to 
stripes as precipitation decreases. However, since the predictions of the amplitude equations break 
down outside a small neighbourhood of the Turing point, it is uncertain whether this stable stripes 
solution will manifest in the full system as a successor to spot patterns as precipitation decreases. 
We expect that as we approach the a = 0 boundary of region A, the interval of stability for the 
spots branch will diminish in size, allowing stable stripes to appear as precipitation decreases. 

In regions B and C of Figure 5, a is negative and 6 -|- 2c and b are positive, which allows a 
stable gaps solution to the amplitude equations in a neighbourhood the upper Turing point. This 
predicts pattern sequences that begin with gap patterns, which is consistent with the standard 
sequence. The stripes solution to the amplitude equations can also be stable in region B, as in 
region A, since c — 6 > 0. Here, gap patterns may transition to stripes as precipitation decreases. 
In region C, c — 6 < 0 prevents the stability of stripe steady states. The analysis therefore provides 
no information in this region about the patterns that may follow gaps as precipitation decreases. 

In regions D and E of Figure 5, a, c — b, and 6 -|- 2c are negative, which means small-amplitude 
steady state solutions cannot be stable near the upper Turing point. Regions D and E differ only 
by the branching direction of the always-unstable stripes solution. The stripes solution branches 
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UpperTuring point 





Figure 5: Summary of upper Turing point amplitude equation calculations over a-f and Dh-f parameter 
spaces, along with schematic bifurcation diagrams. The coefficients of the amplitude equations (4) 
are computed, and the curves a = 0, b — c = 0, c + 26 = 0, and 6 = 0 separate the parameter 
spaces into regions labelled A-E. Qualitatively similar bifurcation structures occur within each 
region. In the white region, no Turing points occur on the uniform vegetated steady state of R02 
and no calculations are performed. 


12 






























































Lower Turing point 
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Figure 6: Summary of lower Turing point amplitude equation calculations over a-f and Dh-f parameter 
spaces. The coefficients of the amplitude equations (4) are computed, and the curves c — b = 0, 
b + 2c = 0, and 6 = 0 separate the parameter space into regions labelled F-I, each of which 
exhibits a qualitatively distinct bifurcation structure. Bifurcation diagrams applicable to regions 
F-I resemble diagrams B-E respectively in Figure 5, with the roles of gaps and spots exchanged 
and the solutions reflected so that supercritical branches bifurcate in the direction of increasing 
precipitation. 


towards the Turing instability interval (i.e. stripes bifurcate supercritically) for region D, since 
6 > 0. The stripes solution branches away from the Turing instability interval (i.e. stripes bifurcate 
subcritically) for region E, since 6 < 0. In both D and E, the gaps solution branches away from 
the Turing instability region. Since there are no small-amplitude steady state solutions stable in 
regions D and E, we cannot directly infer from this analysis what patterned states occur near the 
Turing point. Here, patterned states of the full system likely arise from more strongly nonlinear 
behaviour than the states in regions A-C. 

Eigure 6 summarises the results of bifurcation analysis at the lower Turing point. Over the entire 
a-f and D^-f parameter spaces, a is positive, which is a necessary condition for the stability of spot 
solutions to the amplitude equations. Regions E-H all occur in close proximity to the degenerate 
Turing point curve, while region I fills the majority of the space. Stable solutions to the amplitude 
equations occur only in regions E and G. In region E, the spots solution to the amplitude equations 
is stable near the lower Turing point. The stripes solution is also stable away from the Turing point 
in this region, so that spots may transition to stripes as precipitation increases. In region G, the 
spots solution is stable near the lower Turing point, but the small-amplitude stripes solution can 
never be stable. In regions H and I, solutions to the amplitude equations are never stable near the 
lower Turing point, and differ only in the branching direction of the stripes solution. The stripes 
solution branches towards the Turing instability region for region H, and away for region I. 

At both the upper and lower Turing points, our bifurcation analysis cannot provide direct in¬ 
formation about stable patterned states for a large region of the parameter space, where small- 
amplitude solutions are unstable. The central investigation of this paper is whether the a coef¬ 
ficient, obtained via local analysis, contains information about patterned states near the Turing 
points in these other regions. In regions B and G of Eigure 5, it is expected from the analysis 
that gap patterns are stable near the upper Turing point. We surmise that the same is true in 
regions D and E, where no solutions to the amplitude equations are stable, but a < 0. Similarly, 
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our analysis only shows that spot patterns are stable near the lower Turing point in regions F and 
G of Figure 6 . We surmise that spots patterns will be stable near the lower Turing point in regions 
H and I as well, since a > 0 there. We would then expect to see pattern sequences that begin 
with gaps and end with spots in a scenario of decreasing precipitation over time (i.e. analogs of 
the standard sequence) in the region of parameter space where a < 0 at the upper Turing point 
and a > 0 at the lower Turing point. We find that these conditions are satisfied over nearly all 
of the studied a-f and D^-f parameter spaces, excluding only region A of Figure 5. Region A 
lies adjacent to the Turing degeneracy curve, where the two Turing points approach one another 
and thus the quadratic coefficients at the Turing points approach the same non-zero value. In this 
region, a > 0 at the both Turing points, and we expect pattern sequences that begin and end with 
spots. 

4.2. Numerical simulations 

In numerical simulations, we find that the quadratic coefficient a of the amplitude equations (4) 
signals where the standard sequence occurs in the studied parameter spaces of R02. A summary of 
pattern sequences observed in these numerical simulations is shown in Figure 7. These simulations 
were conducted as described in Section 33.3 and the electronic supplementary material at sets of 
parameter values marked by letters, and pattern sequences were identified by visual inspection. For 
comparison, the region of the parameter space where a > 0 at both the upper and lower Turing 
points (i.e. region A of Figure 5) is shaded. Elsewhere, a < 0 at the upper Turing point and a > 0 
at the lower Turing point. From the results of weakly nonlinear analysis described in Section 44.1, 
we expect pattern sequences beginning with spots (in a scenario of decreasing precipitation) to 
appear in the shaded region. Outside the thin shaded region, we expect to see analogues of the 
standard sequence. As expected, only spot patterned states are observed in numerical simulations 
in the shaded region; in addition, analogues of the standard sequence are primarily observed in 
simulations elsewhere. 

Examples of the numerical patterned states using / = 0.2 and different values of logiQ^D^) are 
also shown in Eigure 7. The simulation output is accompanied by lines which plot the locations 
of the upper and lower Turing points, pu and pc respectively, the transcritical point po, and upper 
and lower pattern stability boundaries, pu+ and pi- respectively (defined in Section 33.3). We 
observed only spot-patterned states in the thin shaded region of Eigure 7. Examples of such states 
are shown in simulation output from / = 0.2 and logiQ^Dh) = 0.6. Near the upper Turing point, 
solutions approximately resemble the spots solution to the amplitude equations (4) shown in 

Eigure 3. The profiles of these spot patterns are roughly sinusoidal about the uniform vegetated 
steady state. An example profile is shown in Eigure 8. At lower values of p, spot patterns remain 
stable. The spacing between spots increases, and the individual spots of vegetation become more 
sharply peaked, quickly decaying to zero away from the centre of a spot. An example of a sharply 
peaked profile is also shown in Eigure 8 . Patterns other than spots are not observed in simulations 
conducted in the shaded region. No notable difference in the qualitative appearance of the spot 
patterns was observed as precipitation increased in discrete steps. 

We primarily observed analogues of the standard “gaps —)■ labyrinths —)■ spots” sequence in the 
unshaded region of Eigure 7, which agrees with our expectations from analysis. Examples of this 
sequence in simulation output for / = 0.2 and logiQ^Dh) ranging from 1.0 - 4.0 are shown in 
Eigure 7. The sets of simulations at logiQ{Dh) = 1.0 and log^o(-^/i) = 1-5 use parameter sets from 
regions B and C of Eigure 5 respectively, where gaps solutions to the amplitude equations are 
expected to be stable near the upper Turing point. As precipitation decreases in the simulations, 
patterns resembling the gaps {H~) solution to the amplitude equations are first observed near the 
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Figure 7: Summary of pattern transitions observed numerical simulations over a-f and D^-f parameter 
spaces in R02, along with representative examples of transitions from numerical simulations at 
/ = 0.2 and log]^Q(I?ft,) = 0.6 - 4.0. Number lines plot the relative locations of the upper and lower 
Turing points (p„ and pe respectively), the transcritical point {po), and upper and lower pattern 
stability boundaries {pu+ and pi-) for the example simulations shown. The parameter values 
corresponding to the example simulations are circled. Though p^ and po are nearly coincident, 
the distance between these points is exaggerated to illustrate that pg > po- 
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Figure 8: Example profiles of individual spot patterns taken from numerical simulations at / = 0.2 and 
logio(4)?i) = 0.6 (sinusoidal), logiQ(£)?,) = 2.0 (sharply peaked), and logio(£)/j) = 4.0 (plateau¬ 
like). The example sinusoidal profile comes from a spot-patterned state near the upper Turing 
point. The sharply peaked and plateau-like profiles come from spot-patterned states well below 
the lower Turing points in their respective simulations. 


upper Turing point. Gaps then transition to well-ordered stripe patterns in both sets of simulations. 
As precipitation decreases further, stripes become disordered before transitioning to spot patterns. 
The sets of simulations at log^gC-^/i) = 2.0 and log^oC-^/i) = 4.0 both use parameter sets from region 
E of Figure 5, where no small-amplitude patterns are stable near the upper Turing point. Still, 
gaps are observed in numerical simulations near the upper Turing point in both sets of simulations. 
As precipitation decreases, gaps transition directly to disordered labyrinthine stripes. These stripes 
eventually transition to spots, which take on non-sinusoidal profiles as shown in Figure 8 . Hysteresis 
occurs in the points of transition between pattern morphologies, and this hysteresis is larger in the 
simulation with the larger value of Dh. The transitions between gaps and labyrinths occur at a 
lower value of p when decreasing precipitation compared to increasing precipitation. The same 
applies to the transition between labyrinths and spots. 

The simulation examples in Figure 7 demonstrate a trend of increasingly nonlinear behaviour as 
Dh increases (see also [31, 32] for similar behaviour in related systems). This trend is generally 
representative of what we observe in the other simulations when parameters are varied away from 
the Turing degeneracy curve, e.g. when a decreases. One aspect of the increasing nonlinearity can 
be accounted for via our bifurcation analysis. The regions in Figure 5 order the parameter spaces 
by nonlinearity at the upper Turing point. Moving away from the Turing degeneracy curve, the 
amplitude equations first predict stable weakly nonlinear patterns in regions A-C, and then imply 
strongly nonlinear patterns in regions D and E since small-amplitude patterns are unstable. This 
manifests in simulations as small-amplitude sinusoidal patterns occurring near the upper Turing 
point when parameters are near to the Turing degeneracy curve, sharply peaked patterns occurring 
beyond this, and plateau-like patterns occurring when very far away from the curve. Other aspects 
of increasing nonlinearity are apparent in certain qualitative behaviours observed in the simulations. 
Patterned states increase in their disorder and begin to exhibit coexistence as distance from the 
degeneracy curve increases. The interval of pattern stability, (p£-,pu+), increases in length as well. 
In particular, decreases to extend the length of the interval {pi-,pi). As this interval increases 
with increasing distance from the Turing degeneracy curve, the transition point between stripe and 
spot states decreases to lower values of precipitation. This causes spot patterns to remain stable 
at values of precipitation well below the lower Turing point. This implies the stabilisation of a 
strongly nonlinear patterned state far from the Turing instability interval. 

In addition to the standard sequence, we observed a few instances of “stripes —>• spots” sequences 
in the unshaded region of Eigure 7. We determined that these are actually instances of the standard 
sequence, where gaps do not appear in simulations. Our bifurcation analysis indicates that gap 
solutions to the amplitude equations are stable only very near to the Turing point in these parameter 
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sets. Because our numerical procedure increments precipitation in discrete steps of fixed size, the 
gaps branch may be bypassed. To test whether gaps can exist stably for parameter sets where 
“stripes —?■ spots” transitions are observed, we conducted additional numerical simulations, which 
are described in the electronic supplementary material. In these simulations, gap patterns were 
assessed to be stable. 

We also observed time-varying spiral wave patterns in one instance of the numerical simulations, 
at / = 0.1 and logiQ{Dh) = 0.5. We ran additional simulations, described in the electronic sup¬ 
plementary material, at nearby parameter values and found that spirals patterns are confined to 
values of that are smaller than typically considered ecologically applicable. To our knowledge, 
there are no previous reports of spiral wave patterns occurring in R02 or any other vegetation 
model. However, we remark that the waves observed in R02 resemble spiral wave patterns ob¬ 
served in other reaction-diffusion contexts such as chemical reaction systems [33] and models of 
phytoplankton dynamics [34]. 


4.3. Quadratic coefficient analysis 

In Appendix A, we derive a closed-form expression for the quadratic coefficient a. This expression 
involves derivatives of nonlinear terms in R02, the infiltration term I{n)h and the growth term 
G{w)n, where /(n) = a{n + f)/{n + 1) and G{w) = w/{w + 1). This expression also involves 
the null vector {Hi, Wi, Ni)'^ and left null vector {Hi, Wi, Ni) of the system linearised about the 
uniform vegetated steady state {ho, wq, hq) (2), both of which are defined in Appendix A. Explicitly, 


a = (^r{no)HiNi + h"{no)hoN^^ (Wi - Hi) 

+ (g'{wo)WiNi + ]^G''{wo)nowi^ (W - 7 W 1 ) . (5) 


The infiltration function I{n) is an increasing concave-down function of n, and so I'{n) > 0 and 
I"{n) < 0. The growth function G{w) is similarly increasing and concave-down with respect to w. 

In a natural limit of the parameters which corresponds to ^ a, we find in Appendix C that a 
negative term dominates quadratic coefficient a (5) at the upper Turing point and a positive term 
dominates at the lower Turing point. These arise through distinct scalings of the critical wavenum¬ 
ber and the onset parameter values with the quantity 5 = I{no)/Dh oc a/Dh at the different 
Turing points, and are calculated in Appendix B. These scalings at onset in turn determine distinct 
scalings for the left and right null vector components. At the lower Turing point, for 5 sufficiently 
small, we find that {Hi,Wi,Ni)'^ = (0(1), 0(1), 0(1))^ and {Hi,Wi,Ni) = {0{5),0{5),0{1)). 
Additionally we find that no = 0{5). Then at the lower Turing point. 


a = + 0(6). 

Ni 


where Qj and W are the positive scaling constants given by (19). Thus a is positive at leading 
order at the lower Turing point. This corresponds to spot patterns. Sufficiently far from the 
degeneracy at which the two Turing points merge, the upper Turing point has the opposite sign. 
Specifkally, at the upper Turing point, we find that {Hi, Wi,Ni)'^ = (0(5^/^), 0(5^/^), 0(1))"^ and 
{Hi, Wi, Ni) = (0(5^/^), 0(1), 0(1)). Then at the upper Turing point 


G'{wo)I''{no)hono ^ 
n + -fG'{wo)no 
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This expression is negative at leading order since I"{n) < 0, and thus corresponds to gap patterns. 

In order for this result to hold, we must also be sufficiently far from a Turing degeneracy. The 
quadratic coefficient at the two Turing points takes on the same sign in a neighbourhood of the 
degeneracy, marked by the shaded region in Figure 7 (see also the discussion of region A in Sec¬ 
tion 44.1). The conditions that ^ a and that we are sufficiently far from a Turing degeneracy 
are not independent. As / increases, the Turing degeneracy occurs at smaller values of 5. This 
can be observed in Figures 5-7, where a decreasing and increasing correspond to smaller 5. At 
smaller S, the regions where a takes the same sign at both Turing points diminish, and less distance 
from the Turing degeneracy is necessary for the scaling results given above to hold. This can be 
seen in the narrowing of the shaded regions in Figure 7 with increasing /. 

5. Discussion 

We find that the quadratic coefficient a from a bifurcation analysis divides the studied parameter 
spaces of the Rietkerk et al. [14] (R02) model into two regions. In a thin region of the parameter 
space adjacent to the degenerate Turing point curve, where the Turing points are very close to each 
other, a is positive at both points. Correspondingly, we observe only spot patterns in numerical 
simulations. Elsewhere a takes opposite signs at the two Turing points. When this happens, we 
primarily observe the standard sequence. This strongly suggests that the a coefficient resulting 
from weakly nonlinear analysis holds predictive value for assessing the nonlinear behaviour of the 
system. Specifically, it appears to serve as a proxy for the sequence of nonlinear patterns that will 
manifest for any parameter set of R02. 

Since a is computed analytically, it is possible to trace the influence of model terms and parameter 
values on the sign of a, and thus on the sequence of patterns that are predicted. This presents an 
approach for comprehensively exploring the full 7-dimensional parameter space of the R02 model. 
Our analysis of the quadratic coefficient in the appendices shows that in any parameter regime 
where the surface water diffusion rate Dh is sufficiently large compared to the infiltration rate a, 
the quadratic coefficient takes values consistent with standard sequence. 

We believe that Dh ^ a is an appropriate limit for the R02 model. We have treated the 
parameters Dh > 1, a > 0 and infiltration feedback strength 0 < / < 1 as essentially unconstrained 
in this study. Our analysis shows, however, that the Turing point calculation is invariant to fixed 
values of the ratio afDh (see Appendix B), and that Turing bifurcations only occur for a/Dh 
sufficiently small. This can be seen in Figures 5-7, where Turing bifurcations only occur for a 
sufficiently small with Dh = 10^, and for Dh sufficiently large with a = 0.4. The amount of 
separation between a and Dh that is required for Turing points to occur depends on /, which can 
also be seen in Figures 5-7. For values of / greater than 0.6, holding other parameters fixed at 
default values, a and Dh must be separated by at least two orders of magnitude for Turing points 
to occur. For any value of / G (0,1), not much additional separation between a and Dh is required 
for the quadratic coefficient a to take opposite signs at the Turing points and for the standard 
sequence to occur. 

Moreover, the limit where Dh ^ a is ecologically relevant. For all soil types, we expect Dh to be 
large, because it represents the ratio of surface water diffusion to biomass dispersal, which occur 
at quite different scales. Additionally, a global study of the factors associated with the existence 
of vegetation patterns by Deblauwe et al. [35] finds that patterns favour environments with non- 
sandy soils, where relatively small rates of infiltration allow for substantial water redistribution 
via surface runoff. For such non-sandy soils, Rietkerk et al. [14] estimate an infiltration rate and 
surface diffusion rate that results in a/Dh ~ 5 x 10“^. Holding the other model parameters fixed 
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at their default values, this corresponds to the standard sequence for 99% of the range of / over 
which patterns are present. 

Although vegetation patterns tend to occur in non-sandy soils, a notable exception is the fairy 
circle phenomenon in Namibia. These patterns occur in a sandy soil environment marked by a high 
infiltration rate [36]. In empirical comparisons between clayey soils and sandy soils, infiltration 
rates differ only by a factor of approximately 10-20 [37]. In the Namibian system specifically, 
estimates for the infiltration rate parameter used in the model by Gilad et al. [4, 11, 22] are at most 
one order of magnitude larger than those in the non-sandy R02 system^ [38]. Given similar surface 
water diffusion rates to R02, this results in the standard sequence over about 95% of the range of 
/ for which patterns are present in the R02 model {a/Dh ~ 5 x 10“^). 

We conclude that the standard sequence prediction appears robust to parameter variation in the 
R02 model. This conclusion is based on evidence that the quadratic coefficient serves as a proxy for 
the standard sequence, along with our finding that the coefficient takes on values consistent with 
the standard sequence in the ecologically relevant region of the parameter space. The methodology 
presented here provides a way forward for assessing the robustness of the standard sequence in 
other models and formulations. Doing so is an important step towards establishing credibility 
for the standard sequence as an early-warning sign, since other model formulations may lead to 
different pattern sequence predictions. The analysis conducted in this paper is specific to the 
formulation of R02. For instance, R02 is distinctive due to the form of its growth function G, 
which depends only on the soil water field {G = G{w) = w/{w + 1)). This leads to a soil water 
value in the uniform vegetated steady state that is constant with respect to precipitation level, and 
only dependent on the mortality parameter (wq = ^/(l — ^)). This feature simplifies the analysis of 
this system. Modifications to the growth function through the addition of dependence on biomass 
(e.g. G = G{w,n)), such as in the models by Gilad et al. [4, 11, 22] and Klausmeier [12, 21] 
alter the structure of the uniform vegetated steady state [39]. Because the value of the quadratic 
coefficient depends on the form of G, we may not find the standard sequence to be similarly robust in 
other systems with different growth functions. Assessing the robustness of the standard sequence 
to alternative (and possibly more general) model formulations is an important and substantial 
undertaking, and could be a direction for future studies. 

Not all PDF vegetation models feature two Turing points on the vegetated steady state. For 
such models our methodology cannot be applied exactly. The Klausmeier Gray-Scott model [12, 21] 
and Simplified Gilad model [4, 11, 22] (in the parameter regime typically studied) have only one 
Turing point on the vegetated steady state. This corresponds, in some sense, to the upper Turing 
point in R02. In the analysis of R02 conducted here, the information used to divide the parameter 
space into regions with different pattern sequences was from the upper Turing point. Because the 
lower Turing point remained positive throughout the parameter space studied, spot patterns were 
always observed in simulations at low levels of precipitation. It may be true that in vegetation 
models with only one Turing point, the sign of the a coefficient at that point is sufficiently informa¬ 
tive to predict where the standard sequence occurs in the parameter space. This would represent 
a remarkable simplification, permitting the efficient analysis of pattern transitions in other PDE 
vegetation models. 
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Appendix 

A. Quadratic coefficient calculation 

In this appendix, we summarise the calculation of the quadratic coefficient, a, of the amplitude 
equations (4) for the Rietkerk et al. [14] model (R02). We do so in a manner that illustrates the 
role played by the nonlinear functions in R02. 

We expand R02 (1) to quadratic order about the uniform vegetated steady state 
which is a function of precipitation p given by (2). We take H = h — /iq; W = w — tco, N = n — uq: 


— \w\=c\w\ + \l2{H,N)--fG2{W,N) I + .... 

dt \ 1 .. I 1 G2{W,N) 


The linear operator C is 

-/(no) + 


C = 


I{nQ) 

0 


0 -/'(no)/io 

-V - '^G'{wQ)no + I\no)ho - jp 

G\wo)no V2 


( 6 ) 


(7) 


the quadratic order terms l 2 {H,N) and G 2 {W,N) are 


hiH, N) = l\no)HN + ^/"(no)/^oiV^ 
G 2 {W,N) = G'{wo)WN+^G"{wo)noW^, 


( 8 ) 

(9) 


and the ellipsis denotes terms of cubic order in {H,W, N). The linear stability of the uniform 
vegetated state to spatially periodic perturbations with wave number q is determined by substituting 
(//, W, into equation (6) linearised about H = W = N = 0. This gives the 

eigenvalue problem where J{q‘^,p) is the Jacobian matrix 




'-I{no)-Dhq^ 0 -I'{no)ho 

/(no) -n - -fG'{wQ)nQ - D^q'^ /'(no)/io - IP 


( 10 ) 


0 


G'{wo)no 


The Jacobian matrix depends explicitly on the wave number q, as well as on the precipitation 
parameter p through the uniform vegetated steady state {ho,wo,no). We consider p to be the 
bifurcation parameter in this analysis. A Turing point occurs at a parameter value p = Pc and 
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wave number = Qc > 0 for which the maximum real part of an eigenvalue is zero, and all other 
modes are damped. Necessary conditions for a Turing point (qciPc) are given by 


Det{J{q‘^,Pc)) = 0, 


dDet{J{q‘^,Pc)) 

dq^ 


<?=9c 


= 0 . 


( 11 ) 

( 12 ) 


We follow a standard procedure [23-25] to obtain the quadratic coefficient a of the amplitude 
equations (4). We write a small amplitude hexagonal (spots/gaps) solution to ( 6 ) near {q'^,Pc) as 

= e {z{h)f{x) + c.c.) j + e^U 2 + 0{e^), (13) 

where e <C 1 , ti = et, 




{Hi, Wi, Ni)'^ is to be determined, and U 2 is a higher-order term which must be bounded. The wave 
vectors lie on a 2D hexagonal lattice, with qi = qd^, 0), q 2 = Q'c( —1/2, \/3/2), qs = — (qi -|- q 2 ). 
Plugging (13) into ( 6 ) gives the 0{e) equation 


0 = J{ql,Pc) 



Since J{q^,Pc) has a zero eigenvalue, we take {Hi,Wi,Ni)'^ to be the associated right null vector 
to satisfy this equation. We choose the convention that + Wf + Ni = 1 and Ni > 0. 

At O(e^), we have 


/ -h{Hi,Ni) \ 

\h{Hi,Ni)-^G 2 {Wi,Ni) {z{ti)f{x) + c.c.f. 

\ G2(Wi,iVi) ) 

The term {z{ti)f{x) -|- c.c.)^ generates modes with wave vectors of magnitude qc- 

{z{ti)f{x) + c.c.f = 2 z^ + c.c. + ... 

These modes result in secular terms in the solution (13). To eliminate these terms, we apply the 
Fredholm alternative theorem to obtain the solvability condition dzjdti = az‘^, where 

a = l 2 {Ni,Hi) (Wi - Hi) + G 2 {Ni, Wi) (w - 7 W 1 ) , (14) 

and {Hi,Wi, Ni) is the left null vector of J{qdPc)- We choose the convention HiHi -|- IFilTi -|- 
NiNi = 2, which eliminates an overall factor of 2 in (14). 


3z 



= CU2 + 
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B. Turing point calculation and scaling 


In this appendix, we derive approximations for the critical wave number and the onset parameter 
value at the upper and lower Turing points, which give scaling relations with respect to the quantity 
6 = I{no)/Dh, where /(n) is the R02 infiltration function, no is the uniform vegetated equilibrium 
biomass, and Dh is the nondimensional surface water diffusion parameter. 

We introduce the abbreviations Jq = I{no), fiQ = I'{no)ho, ho = G\wo)no, where {ho,wo,no) 
is the uniform vegetated steady state given in (2). In the following analysis, we will treat ho as 
the bifurcation parameter, which is justified by the following observations regarding relationships 
between ho, ho, and the precipitation parameter p: 


I. ho is a linear increasing function of p. Specifically, no = {p — Po)/hh where po = np/{l — p), 
and ho is proportional to no by the constant factor G'(wo) = (1 — /x)^ (wo depends only on 
the parameter p). 

II. ho can be expressed as a function of ho, since ho = /'(no)ho = p{no)I' {no) /1 {no) ■ 

III. We will show that Turing points are confined to regions of the parameter space where ho > 7 /x 
(see (17)). 

IV. ho is a monotonic decreasing function of ho when ho > 'yp- This follows from 


dho dno dho 1 dho 

dho dho dno G'{wo) dno 


( 1 -h)^ 


r(n„)ft„ + r(«,)|j2) 


and 


dho 

dno 


d p{no) 


dno I{no) 
1 




p 


I{noy^ 


( 7 /i - I'{no) ho) = - < 0 for ho > 7 /U. 


/(no) ^ /(no) 

Since I{n) is a concave-down function (i.e. I"{no) < 0 ), it follows that dho/dho < 0 . 


Together, these relationships give a one-to-one correspondence between ho, ho, and precipitation p 
in the parameter regime relevant for this analysis. We note that if ho < 'yp at the lower boundary 
of the domain, ho = 0, it can never exceed 7/i, and no Turing points exist. For Turing points to 
exist, it must be the case that ho > 7 /U at ho = 0. This gives a necessary condition for a Turing 
bifurcation in R02: 


ho(0) 


^h(l - /) 
/(I -h) 


> 7h- 


(15) 


This condition explains the asymptotic approach of the degeneracy curves in Figures 5-7 to / = 
V/{v + 7(1 — p)) = 8/9 as a decreases or increases. 

We write the Jacobian matrix (10) as 


J{q^,ho) 


(-I 0 - 

[ 0 “ 


0 —ho \ 

-7ho - n - ho-jpG 

ho -q^ j 


(16) 
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The determinant of this matrix is 


-DhD^q^ - {IqD^ + -fDh{no + / - (^ 0 ( 7^0 + 1 ^) - Dh{ho - 7^)no) - loJfJ-no. 

The zero-eigenvalue condition (11) and the onset condition (12) are necessary for a Turing point 
(9c)™o)) result in the equations 

DwQc + + 7^0 + ^)Qi + ('^( 7^0 + ^) - (K - 71^)^o) 9c + '^7M^o = 0, (17) 

3D^qj + 2 {6D^ + jh^ + i^) q^ + <^( 7^0 +’^) - iK “ 71^)^o = 0> (18) 


where 

^ and Hq = /lo(fiS)- 
Uh 

We treat 5 as a small parameter by assuming 1, noting that Iq is bounded away from zero 

(Iq £ (a/) a))- The only possible negative term in (17) and (18) comes from the factor —(/iq ~ 7 Ai)) 
and so Turing points are confined to regions of the parameter space where ho > 7 /i (as anticipated 
above in observation III). We note that these equations depend on a and only through the 
quantity 5 oc a/Dh. Therefore the Turing point calculation is invariant to fixed ratios of a/D^. 

We proceed by seeking approximate Turing point solutions {q^^ho) to (17) and (18) in the form 
of asymptotic expansions in 5. At 0(1) in S, (17) and (18) are 

Dnjqc + (7^0 + - (hi - lh)noQc = 0 , 

3D^q^ + 2('yh^ + v)ql - (h^Q - 'yn)h^Q = 0 . 

The only real-valued solution in is = 0, which yields (hj) — ')lj)hQ = 0. It can be shown that 
the equation (/iq — 7 /u)hg = 0 has only two unique solutions, hg = 0 and /ig = Thus there 
are two solutions at this order: (i) qj = q^ = H and hg = hg = 0 , which corresponds to the lower 
Turing point, and (ii) = 9 c = 0 and /ig — 7 ^ = 0 where /ig = /ig = 7 /r, which corresponds to the 
upper Turing point. We calculate corrections to these solutions, which yield distinct scalings for 
the wave numbers and onset parameter values with 5 at the different Turing points. 

For the lower Turing point solution, we assume the leading-order correction takes the form 
( 9 |,hg) = (QjS^''-, N£5^^), and we seek values of /3i and ,02 which achieve a balance between terms 
in (17) and (18). For this correction, equations (17) and (18) become 

Qj 

+ ( 5^1 Qj + = 0, 

35^>^^D^Qj + 2 + 52/3i+/32^jY^ + q2 


We observe that 

• DyjQj, DujQj and are always higher order terms than 5‘^^^vQj. 

• —_ '-^n^NiQj must appear at leading order for real solutions in Qj. 

• jg always a higher order term than — 'j^)N£Qj. 
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• must appear at leading order for solutions with Q'j, N^> 0. 

A balance is achieved by Pi = (32 = 1. Thus {qj, Uq) = {Q^S + 0(6“^), N£5 + 0(6“^)), and Qj and Ng 
are solutions to the equations at 


vQl +(v- {ho - 7 fi)Nej Qj + = 0, 

2 uQj + v - {ho~ = 0 , 


(19) 


where ho = ho{N£6) = z^/u(l — /)/(/(! —/u)) + 0{6). It can be shown that a unique physical solution 
to this system (i.e. Qf, N£ > 0) exists when the necessary condition (15) is satisfied. See also 
Dawes et al. [40] for a comparable scaling of the lower Turing point critical wave number with the 
nondimensional water diffusion parameter in the model by von Hardenberg et al. [1]. 

For the upper Turing point solution, we assume the leading-order correction takes the form 
ql = h^o — 7 /U = Invoking a similar argument as for the lower Turing point correction, 

we balance terms of order ^ fjjgj = Pi = 1/2. Thus {q^iho — 7/i) = 

{Ql5^/^ + 0{5),Hu5^l^ + 0{5)), and Ql and are solutions to the leading order equations at 

o{5y. 

{ifio + i')Qt- fioHuQl + 7/^«o = 0, 

2 (7h^ + i^)Ql- h^Hu = 0, 

where ho{ho) = 'y^ + It can be shown that a unique solution to this system, with Ql > 0, 

always exists. 


C. Quadratic coefficient analysis 

In Appendix B, we derived scaling relations between 5 = I{nQ)/Dh, the critical wave numbers qc 
and the Turing point parameter values in terms of ho = G'{wQ)no. In this appendix, we use these 
scaling relations in an analysis of the terms that are important in setting the sign of a when S is 
small. 

The leading order scaling behaviour of the quadratic coefficient a is determined by the scaling 
of the right and left null vectors, (ihi, ITi, A^i)^ and (ihi, ITi, iVi) respectively. We will use the 
relations above to derive the null vector scalings. We recall that a can be written as 

a = (^r{no)HiNi + ^/"(no)/ioiVf) (iTi - #i) 

+ (g\wo)WiNi + ^G”{wo)noWl^ (iVi - yWi) . (21) 

We obtain the following relations between the right null vector components from the first and third 
rows of the Jacobian matrix (10): 

(^Dhql + /(no)) Hi + IpnQhoNi = 0, 

G'(n;o)noITi - qlNi = 0. 


In addition, we recall the right null vector convention Hf + = 1 with Ni > 0. Together, 
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these equations give 


I'{no)hQ 

^ Dhql + I{no) 

Wi = - — - Ni 

G'{wo)no 


N =(l+ gc , I'jnofhl \ 

^ y G'{woynl {Dhq'^ + I{no))‘^ j 

Similarly we obtain the following relations between the left null vector components from the first 
and second columns of (10): 


+ /(no)) - I{no)Wi = 0, 


v + -Dui/c ^G' 


Ni + G'(u;o)noTl^i = 0. 


Together with the left null vector convention HiHi + WiWi + NiNi = 2, we find 

~ ^ I(no) ~ 

^ Dhqc+Hno) 


VTi = 2iVf 


/-/(no)/'(no)/io , + {1 + Dy,)q‘^ 


{Dhqc +I{no)y G'{wo)no 


Ni = 


u + D^qj 

G'{wo)nQ 


+ 7 ITi- 


At the lower Turing point, we found that qj = = Q\b + 0{5‘^) and G'{wQ)nQ = N^b + 0{b^), 

where Qf and are given by (19). The right null vector at the lower Turing point, (22), is thus 
given by 

« r{no)ho ^ 

irA — Aft I rnix\ 


W{ = ^Ni + 0{b), 


^1= l + ^ + 


/'(no)^/io 


Nj (1 + Q 2 ) 2 /(no )2 


+ 0{b). 


Similarly the left null vector components are 


N[ = 2{N[)-^ + 0{b). 


m = 5-—^{Ni)-^ + 0{b^), 


Substituting these expressions in a (21) gives 


a^^ 2QjG'{wo) Ne^^^b) 

I\£ 
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as the leading order behaviour of the quadratic coefficient at the lower Turing point. Since 
G'(u;o), Ne, N( > 0, a( is positive at leading order, corresponding to a prediction of spot pat¬ 
terns near the lower Turing point. 

At the upper Turing point, given that we are sufficiently far from the degeneracy with the lower 
Turing point, we found that q^ = -|- 0{6), where is determined by (20). Given this, 

the right null vector components, (22), evaluated at the upper Turing point are 






QlHno) 


6^/^+ 0(6), 


Ql 

noG'(wo) 

1 + 0 ( 6 ) 


6^/^+ 0(6), 


The left null vector components (23) are 


Hf 


2G'(wo)no 


Qli^ + iQlG'(wo)no 


6G^ + 0(6), 


2G'{wo)no 
V + 7G'(u;o)no 

2 + 0(6^G-^. 


+ 0(6^G^, 


Putting these expressions in a (21) gives 




G'(wo)I''(no)hQnQ 


Since /(no) is a concave-down function, I”(no) < 0, and thus is negative at leading order at the 
upper Turing point. This corresponds to a prediction of gap patterns near the upper Turing point. 
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Supplementary information 

Assessing the robustness of spatial pattern sequences in a dryland vegetation model 

Kama Gowda} Yuxin Chen} Sarah lams^ Mary Silber^ 


1 Numerical simulation procedure 

We numerically solved the Rietkerk et al. [1] (R02) model using the exponential time differencing 
Runge-Kutta 4 (ETDRK4) scheme [2, 3] modified for 2D systems [4], ETDRK4 achieves pseu- 
dospectral accuracy in space and fourth-order accuracy in time. This scheme alleviates issues of 
stiffness often associated with reaction-diffusion systems [4], allowing R02 to be simulated effi¬ 
ciently. Simulations were run using a time step near the scheme’s empirically derived stability 
limit, since our primary concern is with qualitative aspects of patterned solutions. This scheme 
was implemented in MATLAB and C. 

Eor each set of parameter values marked by a letter in Eigure 7, we simulated R02 over a range 
of precipitation values using a procedure described in detail here. Simulations are initialised with 
spatially random initial conditions with fields taking uniformly distributed values in [1,1.5]. The 
simulation begins just below the upper Turing point, at -|-0.95(ptj —pi)- Then the following 

loop is run to identify the upper bound of pattern stability beginning with /c = 0: 

1. R02 is solved using ETDRK4 until either a steady state stop condition or a uniformity stop 
condition is reached. The end state of the simulation, {h,w,n)^ at is saved. 

2. If the steady state stop condition is reached, precipitation is incremented upward by a small 

value Ap, so that p^+^ = -\- Ap. The procedure returns to step 1, using the saved end state 

{h, w, n)^ as the initial condition for a new simulation. If the spatial uniformity stop condition 
is reached, the loop ends. 

The steady state stop condition occurs when either (a) the root mean square difference between 
the current biomass state and the state 400 time units earlier drops below 10“^, or (b) t = 2 x 10^. 
The uniformity stop condition occurs when the root mean square difference between the current 
biomass state and the mean value of that state drops below 10“^. We infer the upper bound of 
pattern stability, pu+ = P'^, to be the value of precipitation where the uniformity condition is 
reached. 

After the first loop terminates at p^ and the first uniformity condition is reached, precipitation 
is decremented by 2Ap so that p'^+^ = p^ — 2Ap. The previous patterned (i.e. non-uniform) end 
state, {h, w, is used as the initial condition for a new simulation. We run the following loop 

to identify pattern morphologies that occur as precipitation decreases, starting with k = N + 1: 

^Department of Engineering Sciences and Applied Mathematics, Northwestern University, Evanston, IL 60208, USA 
^Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA 
^Department of Statistics, The University of Chicago, Chicago, IL 60637, USA 
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1. R02 is solved using ETDRK4 until a stop condition is reached. The end state of the simulation, 
(/i, w, n)^, is saved. 

2. If a steady state stop condition is reached, precipitation is decremented by Ap so that = 

— Ap. The procedure returns to step 1, using the saved end state {h,w,n)^ as the initial 
condition for a new simulation. If the spatial uniformity stop condition is reached, the loop ends. 

When this loop ends at p^~^^ ^ we infer this point to be the lower bound of pattern stability. 

Pi- = , 

We increment precipitation upward in one final loop to identify any potential hysteresis effects. 
Precipitation is first incremented by 2Ap so that = pN+M _|_ 2Ap. The previous patterned 

(non-uniform) end state {h, w, is used as the initial condition for a new simulation. Then 

the following loop is run, starting with k = N + M + 1: 

1. R02 is solved using ETDRK4 until a stop condition is reached. The end state of the simulation, 
{h, w, n)^, is saved. 

2. If a steady state stop condition is reached, precipitation is incremented by Ap so that = 
p^ + Ap. The procedure returns to step 1, using the saved end state as the initial condition for 
a new simulation. If the spatial uniformity stop condition is reached, the loop ends. 

The end result of this procedure is a series of saved end states at a range of precipitation values 
over the interval p G {pi-,pu+). 

The precipitation iteration step size Ap was chosen based on the distance between the upper and 
lower Turing points, Pu — Pi, so that 30-100 end states were saved per parameter set. The value of 
Ap used ranged between 0.0025 and 0.01. Eor most simulations, a time step size of At = 0.4 was 
used, which was near the empirical stability limit of the scheme for most parameter sets. On certain 
parameter sets, values of At as small as 0.01 were needed for numerical stability. We tested for time 
step size errors using a particular parameter set, with At = 0.4, 0.2, and 0.1, and the end states 
were found to be qualitatively indistinguishable. A preliminary set of simulations was run on all 
parameter sets to identify a simulation domain size for each that permitted at least 7 wavelengths 
of patterns. Square L x L domains with L = 400, 800, and 1600 were used. Corresponding N x N 
grid sizes were used, with N = 64, 128, and 256 respectively. 


2 Investigating the “stripes —^ spots” sequences 

Our bifurcation analysis indicates that gap solutions to the amplitude equations should be stable 
near the upper Turing point (region B of Eigure 5). When a is sufficiently small, as it may be for 
the “stripes —)■ spots” observations, the gaps solution is stable for only a relatively small interval 
of the precipitation parameter compared to the stripes branch. This is illustrated schematically 
in bifurcation diagram B in Eigure 7. Because our numerical procedure increments precipitation 
in discrete steps of fixed size, the gaps branch may be bypassed. Also, since gaps are stable only 
very near to the Turing point for these parameter sets, gaps may only be stable at the critical 
wavelength and may not appear in a domain of arbitrary size. 

To test whether gaps can exist stably for parameter sets where “stripes —?■ spots” transitions are 
observed, we seeded numerical simulations at p = pu with a hexagonal lattice gaps pattern at the 
critical wavelength. This initial condition was perturbed by spatially random noise drawn from 
a uniform distribution on the interval [0,0.1]. An aspect ratio of 1 : \/3 and a domain size that 
permits an integer multiple of the critical wavelength were used. In all simulations, gap patterns 
persisted as steady states through t = 1 x 10® and were assessed to be stable in these instances. 
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Supplementary Figure 1: Spiral wave patterns are observed in simulations at / = 0.1 and logio(F*?i) = 0.5. 

An ordered stripes state transitions directly to spirals at p = 0.414. 


3 Investigating the “stripes —> spirals” sequence 

We observed time-varying spiral wave patterns in one instance of the numerical simulations, at 
/ = 0.1 and logiQ{Dh) = 0.5. States from this simulation are shown in Supplementary Figure 1. 
We verified that the observation of spirals is robust to increased spatial and temporal resolution 
in the numerical scheme. We ran additional simulations at nearby parameter values to determine 
whether spirals are confined to a region of the parameter space. Holding / = 0.1 fixed, spirals 
appear in simulations at logig(F>/i) = 0.6 but not log]^o(-^/i) = 0.7. Holding logio(71^/i) = 0.5 fixed, 
spirals appear at / = 0.16 but not at / = 0.18. We conclude that spirals are confined to smaller 
values of than are typically considered ecologically applicable. 
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